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Abstract 

A three species food chain model is studied analytically as well as numerically. 
Integrability of the model is studied using Painleve analysis while chaotic behaviour 
is studied using numerical techniques, such as calculation of Lyapunov exponents, 
plotting the bifurcation diagram and phase plots. We correct and critically comment 
on the wrong results reported recently on this ecological model, in a paper by Rai 
([1995] "Modelling ecological systems", Int. J. Bifurcation and Chaos 5, 537-543). 
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1 Introduction 



Several studies on chaos in continuous as well as discrete ecological models are available in 
the literature [May, 1976; Gilpin, 1979; Schaffer, 1985; Schaffer & Kot, 1985; Hastings & 
Powell, 1991]. Recently Rai [1995] has presented some studies on an ecological model, where 
he has used mainly Painleve analysis for studying the integrability and has given results of 
some numerical calculations also. The calculations given in that paper is completely wrong 
and the claims he is making based on these wrong results are misleading and confusing. 
Almost all results presented in that paper have been published in another paper by Rai et 
al. [1993]. They claim that they could identify the key parameters controlling the dynamics 
of the system, using Painleve property (PP) analysis. They used singular manifold analysis 
of Weiss, Tabor and Carnevale (WTC) [Weiss et al., 1983] which was originally developed 
for testing necessary conditions for Painleve property of partial differential equations and it 
is a direct extension of Ablowitz, Ramani and Segur (ARS) algorithm [Ablowitz et al., 1980] 
for ordinary differential equations. Here we give the correct P-test results for the model. 
Numerical evidence given in that paper is also inconclusive and wrong. Here we present 
results of some extensive numerical calculations and show the irrelevance of their claims. We 
also mention some aspects of the importance of studies on chaos in ecological models. 

2 The Model 

The three species food chain model is given by the following set of ordinary differential 
equations. 

dX AX{1 X BXY 



dt ' " v K' (X + Dx) 
dY -CY + DXY EYZ 



dt (Y + D 2 ) 

dZ 

— = -FZ + GYZ, 
dt 

where A, B, C, D, E, F, G, Di, D 2 and K are model parameters assuming only positive values. 
The model describes a prey population X and two predator populations Y and Z. The 
predator Y preys on X and the predator Z preys on Y. Here Holling type functional 
response [Murdoch & Oaten, 1975] for the predator population is used. Chaos in a similar 
model has been studied earlier [Hastings & Powell, 1991]. 
By the following scaling transformations, 

X _> Kx, Y - *£y, Z - ^-z, t -> j, (2) 



system (0) becomes 



dT y ' (x + a) 



2 



dy yz 

- = -by + cxy-^—^ (3) 

where a = b = c = d = j^t, e = and / = This scaling reduces the number 
of parameters in the system from ten to six. The interesting fact is that, here, K (carrying 
capacity of the environment) and A (growth rate of X) are absent. Since this is a continuous 
system, it can be easily seen that K is just a scaling parameter of the maximum size of 
population X, which can be arbitrarily fixed at some value. By fixing other parameters, K 
or A may be used as a control parameter, which is the case with most of the parameters. 
Parameter E does not have any effect on the dynamics of the system. 

3 Painleve Test 

Painleve test (P-test) is a widely used technique to find the integrability of a dynamical 
system. It is conjectured that if the complex time solutions of a system of equations are 
free from movable singularities other than poles (Painleve property), then it is integrable. 
Several papers and reviews are available on this technique [Ramani et al., 1989; Lakshmanan 
& Sahadevan, 1993]. 

We apply here the P-test as devised by WTC [Weiss et al., 1983] on system ([!]). For testing 
PP in the case of ordinary differential equations, one can use ARS algorithm [Ablowitz et 
al., 1980] which can be considered as a special case of WTC method. Here we try to expand 
the solutions about an arbitrary singular manifold determined by <f>(t) = 0, in power series. 

oo 

x = $>^' +Q 

3=0 

oo 
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where Xj,yj,Zj and are time dependent quantities. If we replace by r = t — to, (to is 
position of the singular point), and consider time independent coefficients, we obtain ARS P- 
test, which is easier to do. But WTC method is useful in finding Backlund transformations, 
Lax Pairs, special solutions, etc. 

Mainly there are three steps in P-test, viz, finding the dominant behaviour, finding the 
resonance values and checking whether there exist arbitrary expansion coefficients at the 
resonance values. Substituting the j — terms of (ffl) in (TJ) and equating dominant terms 
we obtain, a = (3 = —1 and n = — 2. Using this and inserting the full expansion in ([[]) we 
obtain the following recursion relations. 
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On collecting terms of x r ,y r and z r we obtain the matrix equation, 



where 
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and elements of R r are functions of Xj, y iy z iy (i < r) and and their derivatives. When the 
determinant of the matrix M becomes zero, arbitrary coefficients, x r , y r , or z r can exist. 
Those values of r are called resonances. In [Rai, 1995] this step is wrong and hence they 
got wrong values for the resonances. They have only one resonance r = 2, which is wrong 
because in a 3-d system there will be three resonances corresponding to a generic leading order 
behaviour. The actual resonance values are -1 and 1+ ( KD / A } ± V K ® l A +wkd/a+9 ^ j^ egonance 
value -1 corresponds to the arbitrariness of 0. r = 2 is a resonance only when KD/A = 0, 
and resonances are not integers when KD/A > 0. Hence the system is nonintegrable. In the 
case of weak PP we can allow certain rational resonance values also but here it is not valid 
because the dominant powers are integers [Ramani et al., 1982, Joy & Sabir, 1988]. When 
D — 0, the equations for Y and Z are decoupled from that of X, and there is an integer 
resonance r = 2. This special case also does not have PP because there is only one positive 
resonance value; others are negative. 

At the resonance values arbitrary expansion constants should enter in the expansion for 
integrability. But here we can find all the expansion coefficients using the recursion relations. 

For j = : 



Xq = 



yo = 



2(f>t 



and zq = 



20 
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In [Rai, 1995] these solutions are given incorrectly and the expressions they obtained for 
other Xj,yj, and Zj are also wrong. 
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For j — 1 
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Similarly all Xj, yj and Zj for higher values of j can also be determined using the recursion 
relations. For a general solution of the system, we need 3 arbitrary coefficients, which do 
not exist in this case. Thus the nonintegrability of the system is confirmed. 

According to the expressions in [Rai, 1995, Eq. (25)], x\ is determined and their assertion 
that X is indeterminate and hence it is non Painleve type is misleading. They are getting y 2 
and z 2 in terms of x 2 , and claiming that they are arbitrary because x 2 is arbitrary! Moreover 
they are asserting that the system is nonintegrable because these coefficients are arbitrary; 
but for PP we actually need arbitrary coefficients at the resonance values. They claim that 
X is the variable contributing chaotic dynamics to this system which is otherwise completely 
integrable [Rai, 1995, last sentence of the first paragraph in page 540]. Since the system is 
a coupled set of 3 differential equations, only one of the variables cannot be chaotic. When 
we speak about chaos we describe the phase trajectories. It is obvious from basic dynamical 
systems theory that there cannot be chaos in a 2-dimensional system. They have gone to the 
extent of saying that since Xq contains the parameters A and K, the indeterministic nature 
depends on them and they are the key parameters controlling the dynamics of the system; 
and that is their main result! Though the authors did wrong analysis they could give the 
correct answer to the question regarding the integrability of the system in negative. This 
happened to be correct because most of the dynamical systems are nonintegrable. 



4 Numerical Results 

For studying the chaotic behaviour we have plotted the bifurcation diagram and calculated 
the maximal Lyapunov exponents (LE) for various values of the parameters. Fig. 1 gives 
the bifurcation diagram for the parameter values B — 1, C — 1, D — 0.05, E = 1, F = 
I, G — 0.05, D 1 = 10, D 2 = 10, K = 50, when A is varied from 1.0 to 10.0. These are 
the same parameter values which has been used in the paper [Rai, 1995]. To construct 
the bifurcation diagram we integrated the system using the above parameter values and 
after reaching the attractor (discarded large number of initial points), we plotted successive 
maxima (local maxima) of the Z variable, as a function of A. The bifurcation diagram 
indicates a period doubling route to chaos in the system. For higher values of A, sequences 
of periodic windows and chaos repeat. Though the details are different, period doubling is 
observed in each sequence. It is clear from the bifurcation diagram that chaotic behaviour 
starts before A = 4.0 and it is confirmed by the Lyapunov exponent calculation. We have 
checked our calculations with different initial conditions also. 
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We calculated the LE directly from the equations using the same parameter values. For 
calculating the Lyapunov exponents we have to solve the system alongwith the corresponding 
linearised variational system. The method is available in any book on nonlinear dynamics. 
(See for example [Lichtenberg & Lieberman, 1992]). Rai [1995] used Wolf's code for extract- 
ing the LE from a time series, by taking the X variable as the time series, though Wolf et 
al. [1985] lists a FORTRAN code for calculating the Lyapunov spectrum from a system of 
equations also! Methods of calculating LE from a time series are not used when there are 
equations known for describing the system, but it is used when there is only an experimental 
time series available. Finding LE from a time series is not at all reliable and it is always 
approximate. In the calculation of invariants such as LE from time series there are a lot of 
details to be considered such as the number of data points, time delay used for reconstruc- 
tion, the selection of proper embedding dimension, time used for sampling the data, etc. 
There is a large literature available on this topic of time series analysis. (See for example 
a review on this topic by Abarbanel et al. [1993]). In Fig. 2 we give maximal Lyapunov 
exponent of the system as a function of the parameter A from 3.0 to 10.0, keeping all other 
parameters constant. 

Our numerical calculations show that K and A are not the only parameters determining 
the dynamics of the system but most of them have got some relevance. Detailed studies of 
such models in general, will be reported elsewhere. Rai et al. [1993] claim that they did not 
observe chaotic behaviour by varying other parameters. 

In Fig. 3 we have plotted X — Y projection of the attractor for A = 4.0, with all other 
parameters are kept the same as that given above. This is entirely different from the Fig. 2 
of [Rai, 1995]. Here we may note that they have gone wrong somewhere in the calculation, 
because in their Fig. 2 of X — Y plot at A = 4 the X variable reaches values more than the 
maximum it can attain which is equal to the value of K(= 50.); but it goes upto ~ 200. It is 
interesting to note that the maximal LE we got is one order less than that is given in [Rai, 
1995]. We have done the calculations very accurately in double precision for different initial 
conditions and for different variations and verified our results. 

5 Conclusion 

We have done P-test of the model food chain and found it to be nonintegrable. We have 
explored numerically the chaotic behaviour of the system by calculating the maximal Lya- 
punov exponents, plotting the bifurcation diagram and phase plots. We have shown the 
fallacies of the work and irrelevance of the claim that P-test can be used to identify the key 
parameters determining the dynamics of the system by Rai [1995]. In well known models 
like Lorenz, Rossler, etc., also P-test does not give any idea about the key parameters. In 
any system, usually all parameters have some effect on the dynamical behaviour. Usually 
relevant parameters are chosen by their physical importance. In many systems, dimensions 
of the parameter space can be reduced by proper transformations, using symmetries, physical 
considerations, etc. Singularity structure in the complex time plane is very much related to 
the chaotic dynamics of the system, but to my knowledge no paper has appeared in which 
P-test is used to identify directly or indirectly the key parameters controlling the dynamics 
of the system. 
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Selection of biologically realistic parameter values for the numerical simulation of eco- 
logical models is a difficult problem. Biologically relevant parameter values can be found 
by proper identification of the system and its parameters with natural system for which the 
model is applicable. Studies on chaos such as time series analysis, prediction techniques, and 
modelling the system using the nonlinear dynamical data are useful in this regard. With 
biologically realistic parameters, we can have not only chaotic behaviour but also simple 
fixed point behaviour, limit cycles, etc, depending upon the natural system which is studied. 

Many people misunderstand the importance and the meaning of chaos in deterministic 
dynamical systems; they consider the terminology chaos of dynamical systems theory for its 
literary meaning of total disorder or confusion. We can study chaos and even characterize 
it; moreover there is a possibility of checking it with the original natural system. There is a 
kind of order in chaos that is what we are interested in and of course the possibility of short 
term prediction is always there since it is a deterministic system [Schaffer 1985]. Modelling 
of ecological systems and comparing it with actual ecological data help us to understand, 
control, predict, etc., such systems. It helps us to understand how the ecology is going to be 
affected by various external factors also. 
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Figure Captions 

Fig. 1 Bifurcation diagram for the model. Here the local maxima of Z vs A is plotted. 
Other parameters are kept fixed at the values given in the text. 

Fig. 2 Maximal Lyapunov exponent vs model parameter A. Other parameter values are 
the same as that of Fig. 1. 

Fig. 3 Projection of the attractor at A = 4.0 on to the X — Y plane 
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